Microfluidic rheology of soft colloids above and below jamming 
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The rheology near jamming of a suspension of soft colloidal spheres is studied using a custom 
microfluidic rheometer that provides stress versus strain rate over many decades. We find non- 
Newtonian behavior below the jamming concentration and yield stress behavior above it. The data 
may be collapsed onto two branches with critical scaling exponents that agree with expectations 
based on Hertzian contacts and viscous drag. These results support the conclusion that jamming is 
similar to a critical phase transition, but with interaction-dependent exponents. 
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The concept of jamming [1 aims for a unified under- 
standing of dense soft matter such as granular media, 
colloids, foams, emulsions, and glassy liquids. Such sys- 
tems are considered to be jammed if the relaxation time 
is longer than the observation window, so that the parti- 
cles appear stuck in a fixed configuration and the sample 
appears to have a yield stress. At zero temperature and 
zero applied stress, this transition occurs at a critical 
volume fraction of particles, (j)c, and is known as "point 
J" on the jamming phase diagram [5]. The discovery of 
growing dynamical length and time scales on approach 
to point J suggests that it is a critical point '2 - 6 . 

To extend these ideas to driven systems, the shear 
stress a has been studied theoretically as a function of 
both imposed strain rate 7 and of distance |0 — 0c| to the 
critical packing fraction [TlfTU]. The rheology was found 
to collapse onto separate branches above and below cjjc 
when plotted as cr/\(j) ~ (pd^ versus 'y/\4> — 0c T- For 
Hertzian particles in a viscous liquid, conflicting expo- 
nents have been reported as {A = 1.8±0.1, T ~ 2.4±0.1} 
0, {A = 1.5, F = 2.75} [g, and {A = 2, P = 4} (TU]. 
However there has been no experimental test of these 
predictions. In part this is because most experiments are 
restricted to hard spheres, which cannot be packed above 
(j)c ~ 0.64. Also, conventional rheometers can give mis- 
leading results for yield-stress materials due to wall-slip 
and shear banding. Here, we circumvent both issues by 
using soft gel particles and a custom microfluidic rheol- 
ogy technique. We demonstrate critical behavior by the 
collapse of scaled stress versus strain rate, with exponents 
that can be understood in terms of particle-particle in- 
teraction forces in best agreement with Ref. 10^. 

The experimental system is a dense aqueous suspen- 
sion of colloidal N-isopropylacrylamide (NIPA) gel parti- 
cles [TTl [12] , which are swollen with water and hence are 
nearly index and density matched. By changing the tem- 
perature, the amount of swelling and hence the volume 
fraction can be readily controlled. As characterized by 



FIG. 1: (Color onfine) Characterization of NIPA particles, (a) 
Normalized height of a column of particles vs dimensionless 
combination of angular rotation speed u, initial height He, 
and (7 = 9.8 m/s ; colors distinguish difi'erent He- (b) Elas- 
tic modulus of the particles vs temperature, deduced from 
the initial slopes shown as dashed lines in (a), plus empirical 
fit corresponding to -E ~ 1/V'^ where V is particle volume. 
Particle diameters are indicated along the curve. 



dynamic light scattering for dilute samples, the diameter 
is of order 1 ^m and the volume varies with temperature 
as V{T) = (2.93 /im3)[l - T/(39.6 °C)] for our temper- 
ature range, 19 °C< T < 25 °C. The polydispersity is 
about ten percent. Here we study a single sample with 
particle number density of 4.55 x 10^^/m'^, such that the 
empirical fit translates to ^(r) = 1.34 - T/(29.4 °C) for 
our temperature range. 

The mechanical properties of the particles are charac- 
terized by centrifugal compression; see Ref. [13] for full 
details. At rest, the particles sediment to a total pack- 
ing height He corresponding to random close packing of 
spheres at (j)c ~ 0.64. When spun at angular speed uj, 
the packing compresses to a smaller height H, and the 
volume fraction varies with depth such elastic and cen- 
trifugal forces balance everywhere. Data for H/Hc are 
plotted vs (uj'^Hc/g)'^^^ in Fig. [ip, where g — 9.8 m/s^, 
for three different temperatures and for samples with a 
range of He values. The choice of x-axis is both so that 




FIG. 2: (Color online) The experimental setup, (a) The fluid 
is driven through the microfluidic channel by setting the inlet 
and outlet pressures, (b) Video data of the suspension flow 
is taken at mid-channel height, (c) An image taken from a 
video showing the particles. Example velocity profiles are su- 
perimposed for a Newtonian fluid (red) and for NIPA samples 
at = 0.56 (blue) and (j) = 0.64 (green). 



the data collapse for different He, and so that the initial 
decay is linear for Hertzian spheres and can be used to 
deduce the elastic modulus E of the sphere material [13] . 
This holds well; final results for E are plotted vs tem- 
perature in Fig. [T]d, and will used to scale rheology data. 
The sohd black curves in Fig. [T^ represent fits based on 
a constitutive elastic model where the compressive stress 
is Hertzian at small strains and diverges at a finite strain 
[13]. The striking feature is that these fits all asymptote 
to Ha/Hf. « 0.64 at high centrifugal acceleration. Since 
the average volume fraction is inversely proportional to 
height, then Ha/ He = 4'c/4'a- The asymptotic normal- 
ized compression in Fig. [ik therefore gives an asymptotic 
packing fraction oi (j)a = 1- This suggests that the NIPA 
particles deform under compression without deswelling, 
so that packing fractions above <f>c also may be reliably 
computed from the dilute suspension particle size data. 
The steady shear flow rheology is measured with the 
microfluidic device sketched in Fig. [2^. It consists of a 
rectangular PDMS channel, 25 /j,m wide x 100 /im deep 
X L — 2 cm long, fabricated with standard soft lithogra- 
phy [13] and bonded to a glass microscope slide. The fluid 
is forced through the channel using pressurized air and in- 
let/outlet tubing of sufficient diameter that the imposed 
pressure drop AP occurs only along the length L of the 
channel. Force balance therefore allows the shear stress 
at distance y from the center of the channel to be com- 
puted as 17(2/) = {AP/L)y. The corresponding strain rate 
at y is found by numerical differentiation of the velocity 
profile, 7(2/) = dvx{y)/dy [15]. For this, we collect video 
data with a Phantom CMOS camera (1-10,000 fps) con- 
nected to a Zeiss Axi overt 200 microscope with 100 x ob- 
jective focussed at mid-height. Since the channel is tall, 
the observed flow is equivalent to that between parallel 
plates [15]. An objective-coofing collar (Bioptechs) and 
cooling plate above the sample are controlled to about 
0.1 C in order to vary the volume fraction. An example 
video frame in Fig. [2]; displays bead-scale intensity varia- 
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FIG. 3: (Color online) Shear stress versus strain rate for 
several different volume fractions, as labeled. Symbol types 
distinguish runs at different driving pressures. The dashed 
curves represent fits to the Herschel-Bulkley equation. 



tions, so that Particle Image Velocimetry may be imple- 
mented with custom Lab VIEW code. Example velocity 
profiles are superposed on the still image of Fig. [2]:. Alto- 
gether, for a single pressure drop, the cr{y) and 7(y) data 
may thus be combined to give stress vs strain rate shear 
rheology. The dynamic range is typically two decades in 
7, and may be extended by varying the imposed pressure 
drop. 

This rheology concept has been realized previously 
[TBI fT7] . and is related to experiments [TBti2U] where 
the shape of a velocity profile is used to characterize 
shear rheology. Microfiuidics is an ideal platform, since 
the channels are long compared to width so that en- 
trance/exit effects are easily avoided. And owing to the 
small scale, high strain rates may be achieved at low 
Reynolds numbers, so that inertial flow instabilities are 
avoided. Furthermore, the local strain rate is directly 
measured, and hence no problems arise from wall slip or 
shear banding as typically hamper use of conventional 
rheometers for materials with a yield stress. 

Results for stress vs strain rate are collected in Fig. [3] 
Good agreement is found for multiple pressure drops at 
the same volume fraction. This demonstrates the repro- 
ducibility and level of uncertainty in our data; it also 
implies the absence of non-local effects, by contrast with 
Refs. [T7||20]- Note that the data show a clear distinction 
in functional form above and below 4>c- For low (f>, the 
stress tends towards zero at low strain rates. For higher 
(j), the stress extrapolates toward a non-zero yield stress 
ay. We find qualitatively similar results using a conven- 
tional rheometer. To analyze the flow curves, we first fit 
the stress data to the phenomenological Herschel-Bulkley 
form: 

a = ay[l + {JTf]=ay+Kj^, (1) 

where /3 is the shear-thinning exponent, r is a time con- 
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FIG. 4: (Color online) Fitting parameters vs volume frac- 
tion (j): (a) the exponent /3, (b) the dimensionless yield stress 
Oy/E, and (c) the dimensionless timescale tE/tjo; E is the 
particulate material elastic modulus and r/o = 0.01 g/(cm-s) 
is the fluid viscosity. The large symbols are for the main sys- 
tem of particles, as in Figs. 1 and 3; the small symbols in part 
(a) are for particles about 8 times less massive. 



stant, and K is called the consistency. The quality of 
the fits is satisfactory, as shown by the dashed curves 
in Fig. [3] for </> > 0c- The results for /3 displayed in 
Fig. [4^ exhibit no apparent dependence on 0, and have 
average and standard deviation 0.48 ± 0.03. For other 
microgel systems, Ref. [H] discusses yield stress behav- 
ior and Ref. ^ finds /3 = 0.45, while Refs. [UMI fit 
to forms that cross between different limiting viscosities 
at low and high strain rates. The value /3 = 1/2 is pre- 
dicted near jamming for viscously-interacting athermal 
particles [TU] . For simplicity, and so that K has constant 
units, we henceforth fix /3 = 1/2 and repeat the fits. 

The fitting parameters Gy and r are collected in 
Fig. |4|d-c as a function of (j). Both the yield stress and the 
time constant have been rendered dimensionless by ap- 
propriate factors of the elasticity E of the particulate ma- 
terial and the viscosity 770 of the suspending fluid. This 
also serves to eliminate the spurious ^-dependence orig- 
inating from the variation of E with particle swelling. 
While K is always well-defined, Uy and r exist only above 
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FIG. 5: (Color online) Collapse of stress vs strain rate using 
the critical exponents A and P. The dashed lines are fits to 
the Herschel-Bulkley form with /? = 1/2. The first exponent 
uncertainties are statistical; the second are systematic and 
reflect the allowed range of /3. 



jamming and respectively appear to vanish and diverge 
on approach to <j)c. As shown in the main plots, the re- 
sults may be fitted to power-law forms OyjE ^ {(f) — (fjc) 
and TE/r]o - 1/(0 - 0^)^, giving {0^ = 0.633 ± 0.002, 
A = 2.2 ± 0.4} and {0^ = 0.637 ± 0.002, F = 3.8 ± 0.6}. 
The two values for 0c are in agreement and average to 
0.635 ± 0.003, consistent with random-close packing of 
spheres. Fixing 0c to this value, we plot OyjE and tE jri^ 
vs — 0c on log-log axes as insets in Figs. |4]3-c. These 
demonstrate power-law behavior, and give the final re- 
fined scaling exponents as A = 2.1±0.2 and F = 4.1±0.3. 
However we will conservatively take the final statistical 
uncertainties to be twice as large, as given by fits where 
0c fioats. The systematic errors based on the allowed 
range of /3 are 0.1 and 0.4 for A and F, respectively. Note 
that A = /3F holds within uncertainty, which is required 
so that K remains finite and nonzero at 0c and so that 
at high strain rates the stress scales as {7'Iq'J)^E^^~^^ in- 
dependent of 0. Also, the very same exponents are found 
within experimental uncertainty for NIPA particles about 
8 times less massive [T5j. 

Our experimental value of A agrees with that simu- 
lated in Ref. [5], and our full suite of {/3, A, F} values 
are in remarkably good agreement with those predicted 
in Ref. [10] . The observed value of the yield-stress expo- 
nent may be understood physically in terms of the scaling 
of the shear modulus G and the yield strain 7^. For re- 
pulsive particles with interaction energy proportional to 
overlap raised to the power a, numerical simulations find 
G ^ (0— 0c)"~'^/^; this differs from the naive expectation 
a — 2 due to 0-dependent non-affine motion [21 14] . If the 
yield strain scales as Jy ^ (0 — 0c), and if the yield stress 
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[10) . For Hertzian elastic particles, a = 5/2, this predicts 



A = 2 and T = A//3 = 4 as seen here. 

The "distance" (j) — 4>c to jamming thus controls the 
yield stress ay and the time constant r appearing in the 
Herschel-Bulkley form of stress vs strain rate, Eq. (II]) , ac- 
cording to respective scaling exponents A and T. There- 
fore, for volume fractions above (j)c, the shear rheology 
data should all collapse onto a single master curve when 
plotted dimcnsionlessly as a/{E\(j)—(j)c\^) vs 'r]oj/{E\(t) — 
0(;|'"). This construction and the required collapse for 
(f> > 4>c are demonstrated in Fig. [51 A noteworthy feature 
of this plot is that collapse also occurs for all data below 
jamming, for (f) < (j)c, onto a distinct branch. Note that 
the two branches merge close to where the dimensionless 
scaled stress and strain rate are both near 1, which is 
reassuring. The collapse along a second branch need not 
have happened, and serves to emphasize that behavior is 
controlled by distance to point J - just as second order 
phase transitions are controlled by distance to criticality. 
The quality of the collapse is only slightly better using 
/? = 0.48 rather than (3 = 1/2 [12]. As first discussed in 
Ref. [7] , plots such as Fig. [5] represent the measurement 
of universal scaling functions; here, the distinct branches 
are approximately Herchel-Bulkley above (jjc and approx- 
imately power-law below (pc. Recently it was argued that 
collapse may be unique to Hertzian particles [10] . 

Before closing we note that an alternative collapse pro- 
cedure was also explored, using polynomials rather than 
the Hershel-Bulkley form for stress vs strain rate [T5] . 
The best collapse to two distinct branches is found for 
A = 2.2 ± 0.4, F = 5.2 ± 0.8, and 0c = 0.64 ± 0.01. 
This is consistent with the previous analysis, marginally 
so for F, and indicates the extent to which our analysis 
of scaling behavior is model-independent. 

In summary, we have used a custom microfluidic 
rheometer to obtain reliable stress vs strain rate data 
for thermoresponsive NIPA particle suspensions above 
and below jamming, free from shear-banding and wall- 
slip artifacts. Furthermore we have characterized the 
mechanical properties of the particles themselves using 
centrifugal compression, both to account for the change 
in elasticity with temperature and to demonstrate that 
the particles do not deswell when packed above 4'c and 
hence have a well-known volume fraction. Above jam- 
ming, we find that the yield stress scales approximately 
as (0 — 0c)^. For all volume fractions examined, we find 
that the stress increases with strain rate approximately 
as (77")^" where the time scale t grows on approach to 
(j)c as \(f> — (pc\~'^- This represents the first experimental 
measurement of the full set of shear rheology exponents 
on both sides of the jamming transition. The observed 
scaling of stress as power laws of the "distances" |0 — (/)c| 
and 7 to point-J in the jamming phase diagram supports 
the notion that jamming is similar to criticality at phase 
transitions but in a nonequilibrium system. 
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